#include <cmath>
#include "simplex.h"

static int grad3[][3] = {
	{1,1,0}, {-1,1,0}, {1,-1,0}, {-1,-1,0},
	{1,0,1}, {-1,0,1}, {1,0,-1}, {-1,0,-1},
	{0,1,1}, {0,-1,1}, {0,1,-1}, {0,-1,-1},
};

static int p[] = {
	155, 129, 111, 77, 237, 202, 100, 227, 53, 83, 226, 156, 115, 149, 93, 76,
	116, 204, 1, 84, 186, 32, 114, 175, 136, 157, 50, 49, 14, 103, 10, 154,
	236, 192, 167, 169, 159, 108, 119, 242, 231, 47, 144, 200, 42, 44, 110, 198,
	59, 17, 122, 2, 145, 87, 209, 205, 72, 52, 13, 88, 64, 45, 245, 86,
	30, 41, 67, 90, 151, 147, 140, 123, 38, 194, 134, 81, 184, 128, 6, 80,
	255, 71, 20, 21, 172, 208, 65, 39, 235, 79, 254, 220, 54, 12, 180, 247,
	105, 162, 16, 190, 191, 46, 96, 92, 4, 148, 56, 241, 171, 219, 234, 195,
	19, 218, 101, 73, 201, 27, 107, 9, 85, 33, 215, 183, 112, 230, 253, 203,
	118, 246, 31, 187, 207, 206, 224, 165, 48, 70, 251, 176, 109, 233, 26, 181,
	37, 146, 69, 166, 24, 8, 214, 127, 82, 160, 120, 62, 212, 143, 28, 74,
	91, 161, 133, 60, 152, 150, 168, 104, 25, 163, 57, 43, 238, 185, 131, 68,
	117, 252, 174, 229, 124, 138, 63, 78, 225, 34, 216, 248, 164, 217, 211, 173,
	89, 130, 197, 178, 232, 106, 95, 55, 223, 132, 58, 15, 188, 125, 153, 189,
	139, 182, 141, 66, 244, 177, 137, 221, 61, 158, 213, 99, 5, 179, 240, 121,
	228, 29, 22, 18, 113, 11, 51, 126, 35, 135, 222, 142, 94, 170, 249, 210,
	98, 239, 102, 193, 199, 75, 243, 250, 97, 40, 23, 7, 0, 196, 36, 3,
};

static int _floor(double x) {
	return x > 0 ? (int)x : (int)x - 1;
}

static double dot(int g[], double x, double y) {
	return g[0] * x + g[1] * y;
}

double simplex2d(double xin, double yin) {
	double n0, n1, n2;
	const double F2 = .5 * (sqrt(3.0) - 1.0);
	double s = (xin + yin) * F2;
	int i = floor(xin + s);
	int j = floor(yin + s);

	const double G2 = (3.0 - sqrt(3.0))/6.0;
	double t = (i + j) * G2;
	double X0 = i - t,
		   Y0 = j - t;

	double x0 = xin - X0,
		   y0 = yin - Y0;

	int i1, j1;
	if (x0 > y0) {
		i1 = 1; j1 = 0;
	} else {
		i1 = 0; j1 = 1;
	}

	double x1 = x0 - i1 + G2;
	double y1 = y0 - j1 + G2;
	double x2 = x0 - 1.0 + 2.0 * G2;
	double y2 = y0 - 1.0 + 2.0 * G2;

	int ii = i & 255;
	int jj = j & 255;

	int gi0 = p[(ii + p[jj]) % 256] % 12;
	int gi1 = p[(ii + i1 + p[(jj + j1) % 256]) % 256] % 12;
	int gi2 = p[(ii + 1 + p[(jj + 1) % 256]) % 256] % 12;

	double t0 = 0.5 - x0*x0 - y0*y0;
	if (t0 < 0) {
		n0 = 0.0;
	} else {
		t0 *= t0;
		n0 = t0 * t0 * dot(grad3[gi0], x0, y0);
	}

	double t1 = 0.5 - x1*x1 - y1*y1;
	if (t1 < 0) {
		n1 = 0.0;
	} else {
		t1 *= t1;
		n1 = t1 * t1 * dot(grad3[gi1], x1, y1);
	}

	double t2 = 0.5 - x2*x2 - y2*y2;
	if (t2 < 0) {
		n2 = 0.0;
	} else {
		t2 *= t2;
		n2 = t2 * t2 * dot(grad3[gi2], x2, y2);
	}

	return 70.0 * (n0 + n1 + n2);
}
